Probability-Changing Cluster Algorithm for Two-Dimensional XY and Clock Models 
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We extend the newly proposed probability-changing cluster (PCC) Monte Carlo algorithm to the 
study of systems with the vector order parameter. Wolff's idea of the embedded cluster formalism 
is used for assigning clusters. The Kosterlitz-Thouless (KT) transitions for the two-dimensional 
(2D) XY and g-state clock models are studied by using the PCC algorithm. Combined with the 
finite-size scaling analysis based on the KT form of the correlation length, £ oc exp(c/-^/T/Tk;T — 1), 
we determine the KT transition temperature and the decay exponent rj as Tkt = 0.8933(6) and 
n — 0.243(5) for the 2D XY model. We investigate two transitions of the KT type for the 2D g-state 
clock models with q — 6, 8, 12, and for the first time confirm the prediction of r\ — 4/q 2 at Ti, the 
low-temperature critical point between the ordered and XY-like phases, systematically. 

PACS numbers: 75.10.Hk, 64.60. Fr, 05.10.Ln 



I. INTRODUCTION 

The two-dimensional (2D) XY model shows a unique 
hase transition, the Kosterlitz-Thouless (KT) transition 
0. It does not have a true long-range order, but the 
correlation function decays as a power of the distance 
at all the temperatures below the KT transition point. 
Jose et al M studied the effect of the q-fold symmetry- 
breaking fields on the 2D XY model; this is essentially 
the same as treating the q-state clock model, where only 
the discrete values are allowed for the angle 8i of the XY 
spins such that 



9i = 2-Kpi/q with pi = 0,1,2, ■■■,q — 1. 



(1) 



In the limit q — > oo we get the XY model. It was 
shown that the 2D q-state clock model has two phase 
transitions of the KT type at Ti and T 2 (Ti < T 2 ) for 
q > 4. There is an intermediate XY-like phase between 
a low-temperature ordered phase (T < Ti) and a high- 
temperature disordered phase (T > T 2 ). It was predicted 
that the decay critical exponent rj varies from 77 = 1/4 at 
T 2 to r) = 4/q 2 at Ti. 

Most of the above theoretical analyses relied on the 
renormalization group argument, and they are not ex- 
act. There have been extensive numerical studies on the 
2D classical XY model |H-[L(i. In contrast, only a limited 
number of numerical works have been reported on the q- 
state clock model [ pr|Jl^ ] , and the accuracy was not good 
enough especially for the low-temperature phase transi- 
tion. There have been no systematic studies to check the 
prediction of rj(T\) = 4/q 2 . 

In numerical studies, efficient algorithms are impor- 
tant for getting the necessary information. The cluster 
update algorithms of the Monte Carlo simulation [ fj"3||l4| 
are examples of such efforts, and they are useful for over- 
coming the problems of slow dynamics. Recently we pro- 
posed an effective cluster algorithm, which is called the 
probability-changing cluster (PCC) algorithm, of tuning 
the critical point automatically p5j. It is an extension 



of the Swendsen-Wang algorithm [HJ , but we change the 
probability of cluster update (essentially, the tempera- 
ture) depending on the observation whether clusters are 
percolating or not percolating. We showed the effective- 
ness of the PCC algorithm for the Potts models p5|] , de- 
termining the critical point and critical exponents for the 
second-order phase transition precisely with less numer- 
ical efforts. The PCC algorithm was also applied to the 
2D site-diluted Ising model, where the crossover and the 
self-averaging properties were studied [[L6| . The advan- 
tage of using the PCC algorithm in the study of random 
systems is as follows: The sample-dependent T C (L) for 
the finite system with the linear size L is important for 
taking sample average, and the PCC algorithm is suitable 
for getting the sample-dependent T C (L). 

There are a lot of interesting questions about the ex- 
tension of the PCC algorithm, (i) Can the PCC algo- 
rithm be used for the problem of the vector order pa- 
rameter, such as the XY model? (ii) Can it be applied 
to the analysis of the transition of the KT type? (iii) 
Can it work even if the system shows two or more phase 
transitions? 

The purpose of the present paper is to answer these 
questions. We extend the PCC algorithm so as to treat 
systems with the vector order parameter. The rest of the 
paper is organized as follows. In Sec. II, we formulate 
the extension of the PCC algorithm for the vector or- 
der parameter. In Sec. Ill, the KT transition of the 2D 
XY model is studied with the finite-size scaling (FSS) 
analysis based on the KT form of the correlation length. 
In Sec. IV, we study the KT transitions of the 2D clock 
models. We investigate both phase transitions at T\ and 
T 2 for the q = 6, 8, 12 clock models. The summary and 
discussions are given in Sec. V. 



II. PCC ALGORITHM FOR VECTOR ORDER 
PARAMETER 



Our Hamiltonian is given by 



1 



H 



-J Si ■ Sj 



(2) 



where Si is a planar unit vector, (cos#i, sin#i), at site i; 
9i takes the value of [0, 2ir) for the XY model, and the 
value given in Eq. ([!]) for the q-state clock model. The 
summation is taken over the nearest- neighbor pairs 

In order to extend the PCC algorithm to systems with 
the vector order parameter, we use Wolff's idea of the 
embedded cluster formalism juj. We project the vector 
Si onto a randomly chosen unit vector e.\ and another 
unit vector e 2 , perpendicular to e 1; as 



Si = ei cos < 



e.2 sin i 



(3) 



where fa is the angle measured from the axis of the vector 
e.\. Then, the Hamiltonian, Eq. (0), is rewritten as 



n 



UUl) Jl) , r(2) (2) (2) 
j i { ' , ' j 



with positive effective couplings 



J-p = J | COS ( 



(2) 



J I sin i 



(4) 



(5) 



for two sets of Ising variables {ej } and {q 2 •*}. Formally, 
we can restrict ourselves to the region [0, ir/2) for {fa}, 
and we write the partition function as 



,77/2 

Z= / {dfa} 
Jo 



exp(/3 

{ef'=±l} 



{ e < 2) =±i} 



exp(/3£jfef>ef) (6) 



with (3 = 1/ksT. Then, we can use the Kasteleyn- 
Fortuin (KF) cluster representation for the Ising spins 
p7[ . To make the KF cluster, we connect the bonds of 
parallel Ising spins with the probability 



P 



(1,2) 



= 1 - exp(-2/3^' 2) ). 



(7) 



In the PCC algorithm pSf , the cluster representation 
of the Ising model is used in two ways. First, we flip all 
the spins on any KF cluster to one of two states, that is, 
+1 or —1, as in the Swendsen-Wang algorithm |l3[ . Sec- 
ond, we change the KF probability, Eq. (0), depending 
on the observation whether clusters are percolating or 
not. It is based on the fact that the spin-spin correlation 
function G(ri~rj) becomes nonzero for |r< — rj\ — ► oo at 
the same point as the percolation threshold. For the XY 
model in the embedded cluster formalism, the spin-spin 
correlation function is written as 



where (• • •) represent the thermal average. The function 
Q(r i} rj) is equal to 1 (0) if the sites i and j belong to 
the same (different) cluster, and Aij and Bij are some 
constants. Thus, the system is regarded as percolating 
if or Ising spins are percolating. When treating 
the cluster spin update, one may consider Ising spins of 
a single type projected onto a randomly chosen axis, as 
in the original Wolff's proposal fill . However, we should 
consider Ising spins of both types perpendicular to each 
other for checking the percolation. 

The procedure of Monte Carlo spin update is as fol- 
lows: (i) Start from some spin configuration and some 
value of p. (ii) Choose a unit vector e\ randomly, (iii) 
Construct the KF clusters for and using the prob- 
ability, Eq. (^), and check whether the system is perco- 
lating or not. Flip all the spins on any KF cluster to +1 
or —1 for both and Ising spins, (iv) If the system 
is percolating (not percolating) , decrease (increase) fi by 
A/3 (> 0). (v) Go back to the process ii. 

As A/3 becomes small, the distribution of (3 becomes 
a sharp Gaussian distribution around the mean value 
/3 C (L), which depends on the system size L. We approach 
the canonical ensemble in this limit, and the existence 
probability E p , the probability that the system perco- 
lates, becomes 1/2 at f3 c (L). 



III. XY MODEL 

We have made simulations for the classical XY model 
on the square lattice with the system sizes L =8, 16, 32, 
64, 128, 256, and 512. After 20,000 Monte Carlo sweeps 
of determining f3 c {L) with gradually reducing A/3, we 
have made 10,000 Monte Carlo sweeps to take the ther- 
mal average; we have made 100 runs for each size to 
get better statistics and to evaluate the statistical er- 
rors. As for the criterion to determine percolating, we 
have employed the topological rule [ p5p8[ in the present 
study. The topological rule is that some cluster winds 
around the system in at least one of the D directions in 
Z)-dimensional systems. 

Let us start with the size dependence of the transition 
temperature. We use the FSS analysis based on the KT 
form of the correlation length, 



£ oc exp(c/-\/£) 



(9) 



with t= (T — T KT )/T KT . Using the PCC algorithm, we 
locate the temperature Tkt(^) = l/fcs/?c(^) that the ex- 
istence probability E p is 1/2. Then, using the FSS form 
of E p , that is, E v — E p (£/L), we have the relation 



(Si-S 



.ii 
sin < 



I • ± I ( 2 ) (2) 



A ij -eW(r < ,r i )) + ( J By6( 2 )(r i ,r,)) ) (8) 



Tkt(L) = Tkt 



c 2 T KT 
(ln&L) 2 



(10) 



We plot Tkt(L) as a function of I 2 with I = \nbL for 
the best fitted parameters in Fig. [l]. We represent the 



2 



temperature in units of J/ks- The error bars are smaller 
than the size of marks. Our estimate of Tkt is 0.8933(6); 
the number in the parentheses denotes the uncertainty in 
the last digits. We have estimated the uncertainty by the 
X 2 test of the data for 100 samples. This value is consis- 
tent with the estimates of recent studies; 0.89213(10) by 
the Monte Carlo simulation |^|, and 0.894 by the short- 
time dynamics The constant c, in Eq. (|l^), is esti- 
mated as c=1.73(2). 

Let us consider the magnetization (m 2 (L)) at Xkt(L) 
to discuss the critical exponent rj. In Fig. |], we plot 
(?n 2 (L)) as a function of L in logarithmic scale. We ex- 
pect the FSS of the form (m 2 (L)) cx L - '', but there are 
small corrections. The importance of the multiplicative 
logarithmic corrections were pointed out |^,^9|. Using 
the form 

(m 2 (L)) = AL-" (\nb'L)- 2r , (11) 

we obtain 77 = 0.243(5) and r = 0.038(5). We show the 
fitting curve obtained by using Eq. ( |il|) in Fig. ||. This 
value of r\ is a little bit smaller than the theoretical pre- 
diction, 1/4 (=0.25). Our logarithmic-correction expo- 
nent r is compatible with Janke's result r — 0.0560(17) 
for thermodynamic data [[t9| , but different from the the- 
oretical prediction r = —1/16 0|. 



IV. CLOCK MODEL 

Next turn to the g-state clock model. Because of the 
reflection symmetry, we confine ourselves to the case of 
even q. Then, the same procedure can be used as the 
XY model. One thing we should have in mind is that 
the axis of the vector e\ should be chosen from one of q 
directions in Eq. ([[]) or the middle of two of them. We 
plot the high-temperature transition temperature ^(L) 
of the 6-state clock model as a function of l~ 2 in Fig. 0. 
The estimate of T 2 is 0.9008(6), which is more precise 
than the previous estimates; 0.92(1) |l[ and 0.90 [@. 
The plot of (m 2 (L)) at T 2 (L) as a function of L for the 
6-state clock model is given in Fig. [I| The estimate of rj is 
0.243(5) by the analysis of the multiplicative logarithmic 
corrections, Eq. ([TT|), and the exponent r is estimated as 
0.037(5). 

For the second-order transition, the curves of the ex- 
istence probability E p of different sizes cross at T c as 
far as the corrections to FSS are negligible; this is the 
same as the behavior of the Binder ratio [gfj. For the 
KT transition, however, Tkt is not the crossing point 
but the spray-out point. Therefore, T 2 can be searched 
only from the high-temperature side, and Ti only from 
the low-temperature side. The value of E p at T\ is close 
to 1. In principle, we can use the same procedure as the 
study of T 2 ; we may change the setting value of E p , 1/2, 
to a higher one by introducing a biased random walk. 
However, it is difficult to resolve the size dependence for 



lower temperatures. Therefore, we employ a slightly dif- 
ferent approach for the analysis of the phase transition 
at T\. When judging whether clusters are percolating 
or not, we consider another type of clusters. Instead of 
choosing the vector e\ randomly in Eq. (||), we choose 
the vector ei as 

M=\M\ex (12) 

with M = Si, or more generally we may choose e± 
such that 

M = |M| (e x cos0 + e 2 sin0) (13) 

with some fixed angle (f>. With this choice, the existence 
probability for the percolation of only e^ 1 ) (or e^ 2 )) Ising 
spins holds the same FSS property as the total E p . As 

a result, we can control the value of Ep at T\ so as to 
apply the FSS analysis easily with an appropriate <p. 

The low-temperature transition temperature Ti(L) of 
the 6-state clock model obtained by the above modified 
approach is plotted as a function of l~ 2 also in Fig. || As 
the angle <ft in Eq. (|l3|), we have used 7r/3. Our estimate 
of T\ is 0.7014(11), which is more precis e ag ain than the 
previous estimates; 0.68(2) |y} and 0.75 |j. The plot of 
(m 2 {L)) at Ti(L) for the 6-state clock model is also given 
in Fig. [|. The estimate of r\ is 0.113(3) by the analysis of 
the multiplicative logarithmic corrections, Eq. (pi]), and 
the exponent r is estimated as 0.017(4). The previous 
estimates of r\ are 0.100 |ll|] and 0.15 |Lq |. 

We have also made simulations for q = 8 and q = 12. 
The estimates of the transition temperatures T\ , T 2 and 
those of ?7(Ti) and r}(T 2 ) for (j=6, 8, 12, and 00 (the XY 
model) are tabulated in Table [l| The l/g 2 -dependence 
of transition temperatures and exponents are shown in 
Fig. [jj. There, the exact results for q = 4 are also given; 
that is, the Ising sing ularity at T c = [\n(y/2 + l)]" 1 = 
1.1346 with rj = 1/4. The transition temperature T\ be- 
comes smoothly lower with larger q; in the lowest order 
we find that T\ oc l/<7 2 , which is consistent with the the- 
oretical prediction H. The critical exponent rj at T 2 is a 
universal constant, and compatible with the theoretical 
prediction rj = 1/4. The estimates of the critical ex- 
ponent rj at T\ remarkably coincide with the theoretical 
prediction?/ = A/q 2 ; 1/9=0.111 for g=6, 1/16=0.0625 for 
g=8, and 1/36=0.0278 for q=12. This is the first system- 
atic report of confirming the theoretical prediction. 

V. SUMMARY AND DISCUSSIONS 

To summarize, we have extended the PCC algorithm 
|]l5| to the study of the XY and clock models. Wolff's idea 
of the embedded cluster formalism jl4| is used for treat- 
ing the system with the vector order parameter. The KT 
transitions of the 2D XY and clock models are studied 
by using the FSS analysis based on the KT form of the 
correlation length. For dealing with the low-temperature 
transition temperature, T\, we have employed a slightly 



3 



modified algorithm. Investigating the 5 = 6,8, 12 clock 
models, we have systematically confirmed the prediction 
of rj{Ti) = 4/q 2 . We have shown that small logarithmic 
corrections are present in the KT transitions. The sign 
of the logarithmic-correction exponent r is positive for 
all cases of the XY model and the clock models at both 
Ti and T2, which is compatible with Janke's result Jl9| , 
but different from the theoretical prediction r = —1/16 
The present precise numerical results may stimulate 
the refined renormalization-group study of the KT tran- 
sitions. 

In the previous numerical studies of the KT transi- 
tions, one might resort to a bigscale calculation using an 
extensive computer resource [pi, or one might use some 
special boundary conditions B. It is due to the sub- 
tlety of the KT phase transitions; that is, Tkt is not 
the crossing point but the spray-out point of the exis- 
tence probability or the Binder parameter. Moreover, 
the low-temperature transition T\ for the clock model is 
difficult to study because the system is nearly ordered; it 
is more difficult for larger q. We should stress that using 
the present efficient method of numerical simulation, we 
have successfully made a systematic study of the XY and 
clock models with much less efforts. 

Our formalism of the vector order parameter is easily 
extended to the general 0(n) model, where the percola- 
tion of n types of Ising spins will be considered. Then, 
more problems of interest can be studied by the PCC 
algorithm. The PCC algorithm can be also applied to 
the quantum Monte Carlo simulation with the cluster 
algorithm p^,B2[. It will be interesting to compare the 
present result with the quantum XY model. 
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FIG. 1. Plot of Tkt(L) of the 2D XY model for L = 8, 16, 
32, 64, 128, 256, and 512, where I = InbL. 
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FIG. 2. Logarithmic plot of (m 2 {L)) at T K t{L) of the 2D 
XY model for L = 8, 16, 32, 64, 128, 256, and 512. 



FIG. 4. Logarithmic plots of (m 2 (L)) at Ti(L) and T 2 (L) 
of the 2D 6-state clock model for L = 8, 16, 32, 64, 128, 256, 
and 512. 
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FIG. 3. Plots of T,{L) and T 2 (L) of the 2D 6-state clock a function of w a for the 2D tate dock modeL 



model for L = 8, 16, 32, 64, 128, 256, and 512, where Z = InfeL. 



TABLE I. Transition temperatures and exponents r\ for the 
2D q-state clock model. 





T 2 


V(T2) 


Ti 




g = 6 

q = 8 

q = 12 

XY (g = oo) 


0.9008(6) 
0.8936(7) 
0.8937(7) 
0.8933(6) 


0.243(4) 
0.243(4) 
0.246(5) 
0.243(4) 


0.7014(11) 

0.4259(4) 

0.1978(5) 


0.113(3) 

0.0657(2) 

0.0270(5) 
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